Local Linearization - Runge Kutta Methods: a class of A-stable explicit 

integrators for dynamical systems 



H. de la Cruz a < b <*, R.J. Biscay c < d , J.C. Jimenez d , F. Carbonell 

a IMPA, Estrada Dona Castorina 110, Rio de Janeiro, Brasil 
° Universidad de Ciencias Informdticas, La Habana, Cuba 
c CIMFAV-DEUV, Facultad de Ciencias, Universidad de Valparaiso, Chile 
d Instituto de Cibernetica, Matemdtica y Fisica, Calle 15 No. 551, La Habana, Cuba 
c Montreal Neurological Institute, McGill University, Montreal, Canada 



Abstract 

A new approach for the construction of high order A-stable explicit integrators for ordinary differential 
equations (ODEs) is theoretically studied. Basically, the integrators are obtained by splitting, at each time 
step, the solution of the original equation in two parts: the solution of a linear ordinary differential equation 
plus the solution of an auxiliary ODE. The first one is solved by a Local Linearization scheme in such a way 
that A-stability is ensured, while the second one can be approximated by any extant scheme, preferably a 
high order explicit Runge-Kutta scheme. Results on the convergence and dynamical properties of this new 
class of schemes are given, as well as some hints for their efficient numerical implementation. An specific 
scheme of this new class is derived in detail, and its performance is compared with some Matlab codes in 
the integration of a variety of ODEs representing different types of dynamics. 

Keywords: Numerical integrators, A-stability, Local linearization, Runge Kutta methods, Variation of 
constants formula, Hyperbolic stationary points. 
MSC: 65L20: 65L07 



1. Introduction 



It is well known (see, i.e., |ll|, l57[) that conventional numerical schemes such as Runge-Kutta, Adams- 
Bashforth, predictor-corrector and others produce misleading dynamics in the integration of Ordinary Dif- 
ferential Equations (ODEs). Typical difficulties are, for instance, the convergence to spurious steady states, 
changes in the basis of attraction, appearance of spurious bifurcations, etc. The essence of such difficulties 
is that the dynamics of the numerical schemes (viewed as discrete dynamical systems) is far richer than that 
of its continuous counterparts. Contrary to the common belief, drawbacks of this type may not be solved by 
reducing the step-size of the numerical method. Therefore, it is highly desirable the development of numer- 
ical integrators that preserve, as much as possible, the dynamical properties of the underlaying dynamical 
system for all step sizes or relative big ones. In this direction, some modest advances has been archived by 
a number of relative recent integrators of the class of Exponential Methods, which are characterized by the 
explicit use of exponentials to obtain an approximate solution. In fact, their development has been encour- 
aged because their capability of preserving a number of geometric and dynamical features of the ODEs at 
the expense of notably less computational effort than implicit inte grat ors. This have become feasible due 



to advances in the computation of matrix exponentials (see, e.g., [26(, [56|, [ljj, [12], [25]) and multiple 
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integrals involving matrix exponentials (see, e.g., tali J59l | ). S ome instances of this type of integrators are the 
methods known as exponential fitting [47[ , [10| , [60 1, p . 13311 , exponential integrating factor |46|. exponential 
integrators f2" ^.[3~i1j , exponential time differencing 13j ]T [45| . truncated Magnus expansion |36|. [a] , truncated 
Fer expansion [61( (also named exponential of iterated commutators in [35|), exponential Runge-Kutta 



29(, some schemes based on versions of the variation of constants formula (e.g., [50|, [37|, [34 1 



local linearization (see, e.g., [53j], 54 1, 41 1, Q), and high order local linearization methods 



32]. 



The present paper deals with the class of high order local linearization integrators called Local Linearization- 
Runge Kutta (LLRK) methods, which was recently introduced in [l4j as a flexible approach for increasing the 
order of convergence of the Local Linearization (LL) method while retaining its desired dynamical properties. 
Essentially, the LLRK integrators are obtained by splitting, at each time step, the solution of the underlying 
ODE in two parts: the solution v of a linear ODE plus the solution u of an auxiliary ODE. The first one 
is solved by an LL scheme in such a way that the A-stability is ensured, while the second one is integrated 
by any high order explicit Runge-Kutta (RK) scheme. Likewise Implicit-Explicit Runge-Kutta (IMEX RK) 
and conventional splitting methods (see e.g. (48|, [l[), the splitting involved in the LLRK approximations 
is based on the representation of the underlying vector field as the addition of linear and nonlinear com- 
ponents. However, there arc notable differences among these methods: i) Typically, in splitting and IMEX 
methods the vector field decomposition is global instead of local, and it is not based on a first-order Taylor 
expansion, ii) In contrast with IMEX and LLRK approaches, splitting methods construct an approximate 
solution by composition of the flows corresponding to the component vector fields, iii) IMEX RK methods 
are partitioned (more specifically, additive) Runge-Kutta methods that compute a solution y = v + u by 
solving certain ODE for (v, u), setting different RK coefficients for each block. LLRK methods also solve 
a partitioned system for (v,u), but a different one. In this case, one of the blocks is linear and uncoupled, 
which is solved by the LL method. After inserting the (continuous time) LL approximation into the second 
block, this is treated as a non- autonomous ODE, for which any extant RK discretization can be used. On 
the other hand, it is worth noting that the LLRK methods can also be thought of a flexible approach to 
construct new A-stable explicit schemes based on standard explicit RK integrators. In comparison with the 
well known Rosenbrock [j], (HBJ and Exponential Integrators [27j.(29j the A-stability of the LLRK schemes 
is achieved in a different way. Basically, Rosenbrock and Exponential integrators are obtained by inserting 
a stabilization factor (1/(1 — z) or (e z — l)/z, respectively) into the explicit RK formulas, whose coefficients 
must then be determined to fulfil both A-stability and order conditions. In contrast, A-stability of an LLRK 
scheme results from the fact that the component v associated with the linear part of the vector field is 
computed through an A-stable LL scheme. Another major difference is that the RK coefficients involved 
in the LLRK methods are not constrained by any stability condition and they just need satisfy the usual 
order conditions for RK schemes. Thus, the coefficients in the LLRK methods can be just those of any stan- 
dard explicit RK scheme. This makes the LLRK approach greatly flexible and allows for simple numerical 
implementations on the basis of available subroutines for LL and RK methods. 

in [3, [Hi a number of numerical simulations were carried out in order to illustrate the performance 



of the LLRK schemes and to compare them with other numerical integrators. With special emphasis, the 
dynamical properties of the LLRK schemes were considered, as well as, their capability for integrating some 
kinds of stiff ODEs. For these equations, LLRK schemes showed stability similar to that of implicit schemes 
with the same order of convergence, while demanding much lower computational cost. The simulations 
also showed that the LLRK schemes exhibit a much better behavior near stationary hyperbolic points 
and periodic orbits of the continuous systems than others conventional explicit integrators. However, no 
theoretical support to such findings has been published so far. 

The main aim of the present paper is to provide a theoretical study of LLRK integrators. Specifically, the 
following subjects are considered: rate of convergence, linear stability, preservation of the equilibrium points, 
and reproduction of the phase portrait of the underlying dynamical system near hyperbolic stationary points 
and periodic orbits. Furthermore, unlike the majority of the previous papers on exponential integrators, 
this study is carried out not only for the discretizations but also for the numerical schemes that implement 
them in practice. 

The paper is organized as follows. In section 2, the formulations of the LL and LLRK methods are briefly 
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reviewed. Sections 3 and 4 deal with the convergence, linear stability and dynamic properties of LLRK 
discretizations. Section 5 focuses on the preservation of these properties by LLRK numerical schemes. In 
the last section, a new simulation study is presented in order to compare the performance of an specific 
order 4 LLRK scheme and some Matlab codes in a variety of ODEs representing different types of dynamics. 



2. High Order Local Linear discretizations 

Let T> C M. d be an open set. Consider the d-dimensional differential equation 

^ = f(i,x(t)), te[t ,T] (l) 

x(t ) = x , (2) 

where Xo G V is a given initial value, and f : [to,T] x V — > R d is a differentiable function. Lipschitz and 
smoothness conditions on the function f are assumed in order to ensure a unique solution of this equation 
in V. 

In what follows, for h > 0, (t)h will denote a partition to < t% < ... < t/v = T of the time interval [to,? 1 ] 
such that 

sup(h n ) < h < 1, 

n 

where /i n = £ n +i — t n for n = 0, N — 1. 

LocaZ Linear discretization 

Suppose that, for each t n G (tk, y n G 2? is a point close to x(t n ). Consider the first order Taylor 
expansion of the function f around the point (t n ,y n ): 

f (s,u) W f(i n ,yn) + fx(*n,yn)(u - y„) + ft(t„,y„)(s - i n ): 

for s G K and u G £>, where f x , and ft denote the partial derivatives of f with respect to the variables x 
and t, respectively. Adopting this linear approximation of f at each time step, the solution of (HJ-© can be 
locally approximated on each interval [t n ,t n +i) by the solution of the linear ODE 

= A„y(<) + a„ (t) , te[t n ,t n+1 ) (3) 

y (tn) = Yn , (4) 

where A„ = f x (t„, y„) is a constant matrix, a n (t) = f t (t„, y n ) (i - + f (*„, y„) -A„y„ is a linear vector 
function of t. According to the variation of constants formula, such a solution is given by 

t-t n 

y (t) = e A " <*-*»> (y„+ y e- A ""a„(^+ u )d M ). (5) 
o 

Furthermore, by using the identity 

A 

e- A » u du A n = -(e- A " A - I), A > (6) 
and simple rules from the integral calculus, the expression ([S]) can be rewritten as 

y(t) =yn + i>(tn,y n ;t-t n ), (7) 
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where 



t tn 



<t>{t n ,yn;t-t n )= J e A " (t *" u) (A n y n +a„ (t n + u))du 

3 ^(t».y«)(t-*«-«)(f(t n , yn )+f t (t n , yn ) u )d u . (8) 





t-t„ 



In this way, by setting yo = x(io) and iteratively evaluating the expression ([7]) at t n+ \ (for n = 
0, 1, . . . , N — 1) a sequence of points y n +i can be obtained as an approximation to the solution of the 
equation ([1])- <I23) - This is formalized in the following definition. 

Definition 1. (f^J, 0/j For a given time discretization (t) h , the Local Linear discretization for the ODE 
(Ql)- (0) is defined by the recursive expression 

Yn+i = y„ + <p(t n ,y n ; h n ) , (9) 

starting with yo = xo . 

The Local Linear discretization is, by construction, A-stable. Furthermore, under quite general 
conditions, it does not have spurious equilibrium points 4l| and preserves the local stability of the exact 
solution at hyperbolic equilibrium points and periodic orbits 41 1, On the basis of the recursion (O 

(also known as Exponentially fitted Euler, Euler Exponential or piece-wise linearized method) a variety of 
numerical schemes for ODEs has been constructed (see a review in 42|, 16|). These numerical schemes 
essentially differ with respect to the numerical algorithm used to compute ([5]), and so in the dynamical 
properties that they inherit from the LL discretization. A major limitation of such schemes is their low 
order of convergence, namely two. 

2.2. Local Linear - Runge Kutta discretizations 

A modification of the classical LL method can be done in order to improve its order of convergence while 
retaining desirable dynamic properties. To do so, note that the solution of the local linear ODE ©-(j!]) is 
an approximation to the solution of the local nonlinear ODE 

dz(t) „ , . „ , 

— ^ = f f,z i , te[t n ,t n+1 ) 

dt 

z{t n ) = y n , 



which can be rewritten as 

dz(t) 



A n z(t) + a„ (t) + g(t n ,y n ;t,z(t)), t £ [t n ,t n+1 ) 



dt 

z {tn) = y„, 

where g(t n , y n ; t, z (t)) = ((t,z (t)) —A n z(t) — a„ (t), and A„, a n (i) are defined as in the previous subsection. 
From the variation of constants formula, the solution z of this equation can be written as 

z (t) = Yll (t; t n , y„) + r (i; t n , y„) , 

where 

t— t n 

yLL(t;t n ,y n ) = e A - ( - t ~ t ^(y n + [ e^"^ (t n + u) du) (10) 
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is solution of the linear equation ([21)- © and 



t-t„ 



r(f;t„,y„)= / e f ^ t - y ")(*- t "-") g (i„,y„;t„+w,z(i„ + u))dw (11) 



o 

is the remainder term of the LL approximation yi,L to z. Consequently, if r K is an approximation to r of 
order k > 2, then y(t) = yll {t; t ni y n ) + r K (t; t n , y„) should provide a better estimate to z(t) than the LL 
approximation y(t) = yll (t; t n ,y n ) for all t e [t n , t n+ i). This motivates the definition of the following high 
order local linear discretization. 

Definition 2. fsdl] For a given time discretization {t) h , an order 7 Local Linear discretization for the ODE 
(Q]j-(0) is defined by the recursive expression 

Yn+x = Yll (t n + h n ;t n ,y n ) + r K (t n + h n ;t n ,y n ) , (12) 

starting with yo — xo, where r K is an approximation to the remainder term such that ||x(t n ) — y n || = 
0(/i 7 ) with 7 > 2, for all t n £ (t) h . 

Depending on the way in which the remainder term r is approximated, two classes of high order LL 
discretizations have been proposed. In the first one, g is approximated by a polynomial. For instance, by 
means of a truncated Taylor expansion [lij or an Hermite interpolation polynomial [32| , resulting in the so 
called Local Linearization - Taylor schemes and the Linearized Exponential Adams schemes, respectivelly. 
The second one is based on approximating r by means of a standard integrator that solves an auxiliary ODE. 
This is called the Local Linearization- Runge Kutta (LLRK) methods when a Runge-Kutta integrator is used 
for this purpose A computational advantage of the latter class is that it does not require calculation 
of high order derivatives of the vector field f . 

Specifically, the LLRK methods are derived as follows. By taking derivatives with respect to t in ([TTjl . 
it is obtained that r (t;t n ,y n ) satisfies the differential equation 

du (t) 

— -Tj- =<i{tn,y n ;t,u(t)), te[t n ,t n+1 ), (13) 

u(t„)=0, (14) 

with vector field 

q(in) y«; s, = f x (t n ,y n )€ + g {t n , y«; s, y„ + 4> (t n , y„; s - t n ) + £) , 
which can be also written as 

q(in,y n ;s,C) = f(s,y„ + (f>(t n ,y n ;s - t n ) +£)- fx(t n ,y n )(f>(t n ,y n ; s - t n ) 

- ft (tn,y n ) (s - f%) - f (t n ,y n ) i 

where 4> is the vector function §E§ that defines the LL discretization (JSJ). Thus, an approximation r K to r 
can be obtained by solving the ODE (|T3)) - (jl4l) through any conventional numerical integrator. Namely, if 
u Jl+ i = u n +A y " (t n , u„; h n ) is some one-step numerical scheme for this equation, then r K (t n + h n ] t n , y n ) = 
Ay™(t n ,0;h n ). 

In particular, we will focus on the approximation r K obtained by means of an explicit RK scheme of 
order k. Consider an s-stage explicit RK scheme with coefficients c = [a], A = [a-ij], b = [bj] applied to 
the equation (IT3l) - (fT4)) . i.e., the approximation defined by the map 

s 

p (t n , y n ; h n ) = h n ^ &? k i) ( 15 ) 
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where 



This suggests the following definition. 



Definition 3. fflal) An order 7 Local Linear-Runge Kutta (LLRK) discretization is an order 7 Local Linear 
discretization of the form US\) , where the approximation r K to the remainder term Ml)) is defined by the 
Runge Kutta formula \15\) . 



3. Convergence and linear stability 

In order to study the rate of convergence of the LLRK discretizations, three useful lemmas will be stated 
first. 

Lemma 4. Let u n+ i = u n + A y?l (i„, u n ; h n ) be an approximate solution of the auxiliary equation d-?ff)) -i[7 
at t = t n+ i G {t) h given by an order 7 numerical integrator, and y n +i the discretization 



Yn+i = yVi + h n F(t n ,y n ; h n ), 



where 



F(s,£;h) = -{<l>(s,^h)+A^s,0;h)} 

with yo = xo. Then the local truncation error L n+ \ satisfies 

L n+ i = ||x(t n+ i;xo)-x(t n ;xo)-/i n r(t„,x(t n ;xo);/i„)|| < Ci(x. )hZ +1 

for all t n ,t n+ i 6 (t) h - Moreover, if F satisfies the local Lipschitz condition 

\\F{s,fr,h)-F{s,Zi;h)\\<B ( with B e > and 6,6 G e(0 C V, (16) 

where e(£) is a neighborhood of £ for each £ C T>, then for h small enough there exists a positive constant 
C2(xo) depending only on xo such that 

||x(t„+i;x ) -y„+i|| < C 2 (x )/i 7 

for all t n+ i e (t) h . 

Proof. Taking into account that 

x(t n+ i;x ) = y ll ; t n , x(t n ; x )) +r(t n + h n : t n , x(t„; x )) , 

where y^L and r are defined as in (fTO)) and (fTTI) . respectively, it is obtained that 



■ (t n +h n ;t n , x(i„; xo)) - A x(t " ;Xo) (t n , 0; h n ) 



where L n+ i denotes the local truncation error of the discretization under consideration. Since r (t n + h n ; t n , x(i„; Xo)) 
is the exact solution of the equation (|13l) - (fT4")l with y„ = x(i„; Xo) at t n +i and u„ + i = u, l +A x ^ ,i;xo ^ (t n , u„; h n ) 
is the approximate solution of that equation at t n+ \ given by an order 7 numerical integrator, there exists 
a positive constant Ci(xo) such that 

i" {t n + h n ;t n ,x(t n ; x )) - A x ^" ;Xo ^(i n , 0; < Ci(x )^ +1 , 

which provides the stated bound for L n +\- 



On the other hand, since the compact set X = {x (t; xo) : t S [to, T]} is contained in the open set T> C R d , 
there exists e > such that the compact set 

A £ = UeR d : min ||£ - x(i;x )|| < e 

L x(t;xo)GAT 

is contained in V. Since F satisfies the local Lipschitz condition ([T5]). Lemma 2 in [52| (pp- 92) implies the 
existence of a positive constant L such that 

[\F(s,^;h)-F(s^ 1 ;h)\\<L 116-611 (17) 

for all £ Hence, the stated estimate ||x(i n+ i;x ) — y„+i|| < C2(x )/i 7 for the global error 

straightforwardly follows from the Lipschitz condition (TIT)) and Theorem 3.6 in [23| . where C2(xo) is a 
positive contant. Finally, in order to guarantee that y n +i £ A £ for all n = 0, ...,iV — 1, and so that the 
LLRK discretization is well-defined, it is sufficient that < h < d, where S is chosen in such a way that 
C 2 (x )^ < e. □ 



Note that this lemma requires of an order 7 numerical integrator for the auxiliary equation ([13]) - (jT4]) . 
For this, certain conditions on the vector field q of this equation have to be assumed (usually, Lipschitz and 
smoothness conditions). The next two lemmas show that the function </>, and so the vector field q, inherits 
such conditions from the vector field f . 

Lemma 5. Let ip(.; h) — j^<f> (.; h). Suppose that 

f e C p+hq+1 ([t ,T] x V,R d ) , 
where p,q € N. Then <p G C p > q > r {[t a , T] xVx M + ,M d ) for all r G N. 
Proof. Let $j be the analytical function recursively defined by 

(^ j (z)-l/j\)/z forj = 1,2,.. 



(z) 

for zgC. Since 



i = o 



for all s G M+ and M G R d x M d (see for instance [Hg), the function 9? can be written as 

<p (r, £; <J) = a (<5f x (r, 0) f (r, £) + tf 2 (<5f x (r, 0) f t (r, «5 

for all r G M, £ G R d and <5 > 0. Thus, from the analyticity of dj and the continuity of f the proof is 
completed. □ 

Lemma 6. Let f and q be the vector fields of the ODEs (Ip-flP and 113\) - [T4\ ), respectively. 

i) There exists e > such that the compact set 



A e = \zeR d : min llx (t) - zll < e 1 
1 te[t D ,T] " J 



is contained in P. Moreover, there exists a compact set /C e included into an open neiborhood of and 
a S e > 0, such that 

x(t) + </>(£, x(f );£)+£ G A, 
for all <5 G [0, 5 e ], £ e /C E and t G [t 0) T]. 
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ii) If f and its first partial derivatives are bounded on [to,T] x T>, and f(t, .) is a locally Lipschitz function 
on T> with Lipschitz constant independent of t, then there exists a positive constant P such that 

||q(t, x(i); t + 6, &) - q(t, x(f); t + <S, <P\\&- 611 

for all 5 G [0, 4], 6,6 € £ e and i G [t , T). 
m) If f e C p {[t ,T] x P,E d ) for somcpe N, then q(t,x(i); •) € C p ([i,i + <5 e ] x K, £ ,R d ) for all £ € [t ,T]. 

Proof. The first part of assertion i) follows from the fact that X = {x (i) : t £ [to,T]} is a compact set 
contained into the open set T>, whereas its second part results from the continuity of <h on [to, T] x A £ x [0, 8 e ] 
stated by the Lemma[SJ Assertion ii) is a straighforward consecuence of Lemma 2 in [52j (pp. 92). Assertion 
Hi ) follows from the definition of the vector field q and Lemma [5j □ 

The next theorem characterizes the convergence rate of LLRK discretizations. For this purpose, for all 
t n E (t) h , denote by 

y«+i = y n + h n Lp 1 (t ni y n ; h n ) (18) 
the LL discretization defined in (fT2"|) . taking r K as an order 7 RK scheme of the form (fTS)) . That is, 



y 7 (in,yn; h n ) — — {(j) (t n ,y n ; hn) + p(t n ,y n ; h n )} , 



where cj> is denned by ©. 
Theorem 7. Suppose that 

Then 



f e C + \[t ,T] x V,R d ). 



(19) 



||x(£„ + h;xo) - x(t„;x ) - hcp y (t n> x(i„; x ); h)\\ < d(x )/i 7+1 , 
and the LLRK discretization U8\) satisfies 

||x(i n+ i;Xo) -y n +i|| < C 2 (xo)/i 7 , 
for all t n ,t n+ i 6 (t) h , where Ci(xo) and C2(xo) are positive constants depending only on Xo. 

Proof. By Theorem 3.1 in [23}], the local truncation error of the order 7 explicit RK scheme ([T5|) for the 
equation (Tl3|) - (ll"4l) with y„ = x(i„;x ) is 



\\u(t n + h)-p(t n ,x(t n ;x );h)\\ <C(x ) h~< +1 , 



(20) 



where 



with 



C(x ) 



1 



(7 + l)!ee[o,i] 



d^ +1 



dri- 



T u{t n + eh) 



i 8 

h 2> 



06 [0,1] 



—k t (6h) 

dp n ; 



k;(0ft) = q ( t n , x(£ n ; xq), t n +Cj9h,9h cufa^h) 

By taking into account that the solution r of (IT3"1) - (|14I) is the remainder term of the LL approximation 
and by setting y n = x(i„;xo) in (|13[) . it follows that 



u (t„ + 6h) = x (i n + x ) - x (t n ; x ) - (t„, x (£„; x ) ; 0/i) , 



(21) 



and so 



dP+ 



T u(t n + Oh) 



d 1 

— q(t„ x(i n ; xq); i„ + 0/i, u(i„ + 9h)) 



where the derivative in the right term of the last expression is with respect to the last two arguments of the 
function q. Condition (fT9"|) . assertion Hi) of Lemma [5] and expression (|2ip imply that q(.,x(.,xo); -,u(.)) G 
C 7 ([io, T), R d ). Hence, there exists a constant M such that 



max 
0e[o,i], t„e[t ,T] 



d'< +1 

^u(t n + 6h) 



< M. 



Likewise, condition fj 19|) and Lemma [6] imply that 



max 

9e[0,l], t„G[t ,T] 



d 7 

-T-ki(fl/i) 

dt~t v ; 



< M. 



Therefore, C(x ) in (j2"0)) is bounded as a function of x G V. 

In addition, Lemma[5]and Lemma 3.5 in [23] combined with assertion Hi) of Lemma [6] imply that cf> and 
p satisfy the local Lipschitz condition (|16p . and so does the function 



<py(tn,y n ;h) = -{(t)(t n ,y„]h) + p{t n ,y„;h)} 



as well. This and Lemma 2] complete the proof. 



□ 



Note that the Lipschitz and smoothness conditions in Lemma|4]and Theorem[7]are the usual ones required 
to derive the convergence of numerical integrators (see, e.g., Theorems 3.1 and 3.6 in [23]). These conditions 
directly imply that smoothness of the solution of the ODE in a bounded domain (see, e.g., Theorem 1 pp. 79 
and Remark 1 pp. 83 in [52jj ) . In this way, to ensure the convergence of the LLRK integrators, the involved 
RK coefficients are not constrained by any stability condition and they just need to satisfy the usual order 
conditions for RK schemes. This is a major difference with the Rosenbrock and Exponential Integrators and 
makes the LLRK methods more flexible and simple. Further note that, like these integrators, the LLRK are 
trivially A-stable. 



4. Steady states 

In this section the relation between the steady states of an autonomous equation 

*^ = f(x(t)), te[t ,T], 
x(to) = x G M d , 



(22) 
(23) 



and those of their LLRK discretizations is considered. For the sake of simplicity, a uniform time partition 
h n = h is adopted. 

It will be convenient to rewrite the order 7 LLRK discretization in the form 



where 



with 



y«+i = y n + hip^yn, h), 

s 



i=l 



(24) 
(25) 

(26) 



ki (£,<5) = q(£; <h6, S a l0 kj (£,5)) 



and 



q (6 s, u) = f (e + <5$(e, 5)f (o + u) - uosm mo - f (o • 



For later reference, the following Lemma states some useful properties of the functions ip-y on neighbor- 
hoods of invariant sets of ODEs. 



Lemma 8. Let S C M. d be an invariant set for the flow of the equation 



Let K, and be, respectively, 



compact and bounded open sets such that E C K, C fl. Suppose that the solution x of 
condition 

x(t; Xo) C f2 for all initial point xo G K and t 6 [fojT], 
and the vector field f satisfies the continuity condition 

f e C 7+1 (ft,K d ). 

Further, let 

yn+i = y« + hip y (y n , h) 
be the order 7 LLRK discretization defined by {21$. Then 



fulfils the 
(27) 

(28) 



z) ip~f — > f and dif^/dy n — > f x as ft — > uniformly in K, 
ii) II (x(to + x o) — Xo)//i — tp 7 (xo, ft) II = 0(ft 7 ) uniformly for xq G /C. 

Proof. According to Lemma 5 in [41|, f e C 7+1 (f2) implies that $/ -4 f and <9($/)/<9£ — > f x as ft — > 
uniformly in /C. 

On the other hand, kj(£, 0) = 0, for all £ £ il and i = 1, . . . , s. Besides, since 



9£ 



3=1 



where 



||(£; <5, u) = f x (£ + «S)f (£) + u) A (£ + 5 )f (f) + u ) 

- <^ (fx(o^(e, *)f (0) - fx (o + fx(e + , *)f (o + u) ^ 



with 



u = S o-ijhj (£, <5) and 
3=1 



du 



3=1 



_5 



9kj(£, 0)/<9£ = for alH = 1, . . . , s. Thus, since each and <9k;/<9£ are continuous functions on £1 x [0, 1], 
it holds that k; -> and 9kj/9£ — > as ft — > uniformly in the compact set /C. Thus, assertion i) holds. 
From Theorem [7] we have 



(x(t + ft;x ) -x )/ft-^ 7 (x ,ft)|| < C(x )ft 7 , 



where 



C(x ) 



1 



max 



(7 + l)!ee[o,i] 
is a positive constant depending of Xq, 



d 7 + 



T u(t + Oh) 



1 - 

H r > l&il max 

-v! ^— ' ee[o,i] 



i = l 



d 7 

— k 4 (#ft) 



i-1 



ki(6>ft) = q ( x(t ; x ); CjOhfih aijkj(0h) J , 

3=1 



i=l,.. 



10 



and u (t + Oh) = x (f + 6h; x ) - x (to', x ) - (to- x (t ; x ) ; Oh). 
Clearly, 



df+ 



T u(s) 



<ri+ 



— pj-(x(s;xo) -x(t ;x ) - <t> (to, x (t ; x ) ; s - t )) 



< 



dt 7 



f(x(s;x )) 



d 7+ 



r </)(to,x(to;xo);s-t ) 



for all s £ [to,^o + ft]- Since x(t;xo) G f2 for all t € [to,? 1 ] and xq e /C C there exists a compact set Ah 
depending of h such that /C C Ah C ^ and x(s;x ) G ^ for all s S [to, to + ft] and Xo € 1C. In addition, 
since condition f g C 7+1 (f2, R d ) implies that there exists a constant M such that 



sup 



d 7 



it is obtained that 



max 
ee[o,i], x £C 



d 7 
di 7 



f(x (t o + 0ft;x o )) 



< M, 



< sup 



d 7 
dt 7 



f(0 



< M. 



Taking into account that <f> and are functions of f , we can similarly proceed to find a bound B > 

r<?!> (to, x (to; Xo) ; s — t ) and || Jj^kj(0ft)||. Hence, we conclude that C(xo) is 

□ 



independent of 6*, Xo for 

bounded on K. by a constant independent of Xq, and so ii) follows. 



dt~< + 1 



4-1. Fixed points and linearization preserving 

Theorem 9. Suppose that the vector field f of the equation i2S\) and its derivatives up to order 7 are 
defined and bounded on R d . Then, all equilibrium points of the given ODE (2I2\) are fixed points of any 
LLRK discretization. 

Proof. Let ip 7 , $ and k^ be the functions defined in expression (|25l) . If £ is an equilibrium point of (|22p. 
then f(£) = and so $(£, ft)f(£) = and kj(£, ft) = for all ft and i = 1, . . . , s. Thus, <p 7 (f, ft) = for all 
ft, which implies that £ is a fixed point of the LLRK discretization □ 

A numerical integrator u n +i = u„ + A (f„, u n ; ft„) is linearization preserving at an equilibrium point £ 
of the ODE (|22p if from the Taylor series expansion of A (t n , •; h n ) around £ it is obtained that 

u„+i - £ = e M ^(u n -0 + 0(K - e|| 2 ). 

Furthermore, an integrator is said to be linearization preserving if it is linearization preserving at all equi- 
librium points of the ODE (i^ . 

This property ensures that the integrator correctly captures all eigenvalues of the linearized system at 
every equilibrium point of an ODE, which guarantees the exact preservation (in type and parameters) of a 
number of local bifurcations of the underlying equation [49[ . Certainly, this results in a correct reproduction 
of the local dynamics before, during and after a bifurcation anywhere in the phase space by the numerical 
integrator. 

In (Zij the linearization preserving property of the LL discretization ^ was demonstrated. This property 
is also inherited by LLRK discretizations as it is shown by the next theorem. 

Theorem 10. Let the vector field f of the equation and its derivatives up to order 2 be functions defined 



id bounded 



Then, LLRK discretizations are linearization preserving. 



Proof. Let £ be an arbitrary equilibrium point of the ODE 
neighborhood of £. 

Let us consider the Taylor expansion of f around £ 



and let the initial condition y„ be in the 



f(y») = fx(0(y»-O + 0(lly»-£ll ) 
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and the LL discretization 

y n +i = y n + h$(y n ,h)f(y n ), 

where <I> defined as in (|2T))) is, according to assertion i) of Lemma 1 in [4l|, a Lipschitz function. By com- 
bining this Taylor expansion with both, the identity © and the Lipschitz inequality ||$(y n) h>) — h)\\ < 
^ ||y« — £|| it is obtained 

M>(£, fc)f (y„) = (e ht *^ I)(y„ - + 0(||y„ - £|| 2 ) (29) 

and 

||(*(y»,/i)-*K,/i))f(yn)|| <^||y„-c|| 2 , (30) 

respectively, where C is a positive constant. 
Now, consider the LLRK discretization 

y«+i = y« + hip^yn, h), 

with (p-y defined as in (|25l) . From the Taylor formula with Lagrange remainder it is obtained that 

||q(y„; h, u)|| - ||f(y„ + $(y„, h)f(y n )h + u) - f x (y„)$(y„, /i)f (yn)*» - f (y„)|| 
< M ||$(y„, h)i(y n )h + uf + ||f x (y„)|| ||u|| , 

where the positive constant M is a bound for ||f xx || on a compact subset Kcl'' such that y„, y n +i,£ £ 
By using ^ and (O it follows that 

||q(y„; h, u)|| < 2M ||u|| 2 + ||f x (y n )|| ||u|| + 0(||y„ - e|| 2 ). 

From the last inequality and taking into account that ki = 0, it is obtained that ||k.2|| < 0(||y n — £|| 2 ). 
Furthermore, by induction, it is obtained that ||kj|| < 0(||y„ — £|| ) f° r a ^ ^ = 1,2, s. From this, (12T)1) and 
(1251) it follows that 

h<p 7 (y n , h) = ( e ^C« - I)(y„ - £) + 0([[y„ - £|j 2 ), 
which implies that the LLRK discretization is linearization preserving. □ 

The next two subsections deal with a more precise analysis of the dynamical behavior of the LLRK 
discretizations in the neighborhood of some steady states. 

4-2. Phase portrait near equilibrium points 

Let be a hyperbolic equilibrium point of the equation (|2"2l . Let X s , X u C R d be the stable and 
unstable subspaces of the linear vector field f x (0) such that M. d — X s © X u , (x s ,x„ ) = x G R d and || x = 
max{||x s || , ||x„||}. It is well-known that the local stable and unstable manifolds at may be represented 
as M s = {(x s ,p(x s )) : x s G fC e ,s} and M u — {(q(x u ), x u )) : x u G ?C e ,«}, respectively, where the functions 
p : K. e>s =K £ flI s ^ /C e!tl = K £ n X M and q : fC e:ll — > £ EjS are as smooth as f , and fC e = {x E R d : ||x|| < e} 
for e > 0. 

Theorem 11. Suppose that the conditions \&7ty -(ES \) of Lemma \^ hold on a neighborhood f2 of 0. TTien 
i/iere ea;?si constants C, e, Eq, ho > suc/i i/iai f/ie /ocaZ stable and unstable manifolds of the order 
7 LLRK discretization \2J$ at are of the form 

M h s = {(x s ,/(x s )) : x s G /C E , S } and M£ = {(^(x u ), x u ) : x u G /C E ,»}, 

where p h = p + 0(/i 7 ) uniformly in K e ,s , and q' 1 = q + 0(/i 7 ) uniformly in IC e . u . Moreover, for any xq G /C e 
and h < ho, there exists zo = zo(xo, h) G K, eo satisfying 

sup{\\x(t n ;x ) -y„(z )|| : x(t;x ) G /C E /or t G [*o,*n]} < C/i 7 . (31) 

Correspondingly, for any zq G IC e and h < ho, there exists xq = xo(zo,/i) G IC £o that fulfils \31)) . where the 
supremum is taken over all n satisfying y 3 (zo) G K, e , J = 0, ... ,71. 
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Proof. Since Q is a neighborhood of the invariant set 0, there exists a constant e > and a compact set 
JC e — {£ 6 K" : ||£|| < e} C f2 such that Lemma [5] holds with /C = IC e . Furthermore, by assertion i) of 
Theorem[5J f (£) = implies <^ 7 (£, h) = for all /i. Thus, the hypotheses of Theorem 3.1 in Q hold for the 
LLRK discretizations, which completes the prove. □ 

Theorem [TT] shows that the phase portrait of a continuous dynamical system near a hyperbolic equi- 
librium point is correctly reproduced by LLRK discretizations for sufficiently small step-sizes. It states 
that any trajectory of the dynamical system can be correctly approximated by a trajectory of the LLRK 
discretization if the discrete initial value is conveniently adjusted. It also affirms that any trajectory of a 
LLRK discretization approximates some trajectory of the continuous system with a suitably selection of 
the starting point. In both cases, these results are valid for sufficiently small step-sizes and as long as the 
trajectories stay within some neighborhood of the equilibrium point. Moreover, the theorem ensures that 
the local stable and unstable manifolds of a LLRK discretization at the equilibrium point converge to those 
of the continuous system as the step-size goes to zero. 

4-3. Phase portraits near periodic orbits 

Suppose that the equation (|22[) has a hyperbolic closed orbit F = {x(i) : t 6 [0,T]} of period T in an 
open bounded set fl C R rf . Let fl be the closure of f2. 

Theorem 12. Let the assumptions l[27\ )- [28\) of Lemma\8\ hold on a neighborhood ofQ. Then there exist 
ho > and an open neighborhood UofT such that the order 7 LLRK discretization 

Yn+i = y n + hip y (y n , h) 

has an invariant closed curve T/,. C U for all h < ho- More precisely, there exist T —periodic functions 
y h : R — > U and — 1 : R — » M for h < ho, which are uniformly Lipschitz and satisfy 

y H (t) + fup 7 (y h (t),h) =y h (a h (t)), t € R 

and 

a h (t) =t + h+ 0{h 1+1 ) uniformly for t G R. 
Furthermore, the curve T% = {y^(^) :tE [0,T]} converges to T in the Lipschitz norm. In particular, 

max ||5E(t) -y h (t) \\ = O(^) 

and 

\\(x-y h )( tl )-(x-y h )(t 2 )\\ 
sup j j > U as h — > U. 

Proof. Since Lemma[5]holds on a neighborhood of f2, it also holds on ft. In addition, Lemmas [5] and [5] imply 
that V'x.G C 2 (f2 x [0, ho]), so d<p 1 /dy n is Lipschitz on ft uniformly in h. Thus, the hypotheses of Theorem 
2.1 in [3j hold for the LLRK discretizations of order 7 > 2, which completes the proof. □ 

Theorem [T^] affirms that, for h sufficiently small, the LLRK discretizations have a closed invariant curve 
Th , i.e., (1 + htp(.; /i))(r/j) = F^ , which converges to the periodic orbit T of the continuous system. 

The next theorem deals with the behavior of the discrete trajectories of LLRK discretizations near the 
invariant curve Th when the ODE ([2^|) has a stable periodic orbit T. For Xo in a neighborhood of F, the 
notations 

W h (x ) = {y„(x ) : n > 0} and w(x ) = {x(t; x ) : t > 0} 
will be used. In addition, 

d(A,B) = maxjsup dist(z, B), sup dist(z,A)} 

z£A zS-B 

will denote the Hausdorff distance between two sets A and B. 
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Theorem 13. Let T be a stable closed orbit of the equation \22jl . Then, under the assumptions of Theorem 
WA there exist ho, a, (3, C and p > such that for h < ho and dist(jXQ,Th) < p the following holds: 

dist(y n (xo),T h ) < C exp(-crf„) dist(x ,T h ) 

and 

dist(y n (x ),w(x )) < C(/i 7 + min{h' y exp(/3f n ), exp(-ai n )}) 
for n > 0. Moreover, for any S > there exist p{5), h{5) > such that 

sup{dist(y n (x ), tu(xo))} < C/i 7 " 5 

n>0 

for h < and dist(xo,r; l ) < Finally, 

d(Wh(x ),w(x a )) ->• as /i ->• 

uniformly for disi(xo, T) < p. 

Proof. It can be proved in a similar way as Theorem 1121 but using Theorem 3.2 in Q instead of Theorem 
2.1. □ 

This theorem states the stability of the invariant curve F/j and the convergence of the trajectories of a 
LLRK discretization to the continuous trajectories of the underlying ODE when the discretization starts at 
a point close enough to the stable periodic orbit V. 

5. A-stable explicit LLRK schemes 

This section deals with practical issues of the LLRK methods, that is, with the so called Local Lineariza- 
tion - Runge Kutta schemes. 

Roughly speaking, every numerical implementation of a LLRK discretization will be called LLRK scheme. 
More precisely, they are defined as follows. 

Definition 14. For an order 7 LLRK discretization 

y n +i = y n + h n ip 7 (t n ,y n ; h n ), (32) 
as defined in 118\) . any recursion of the form 

y„+i = y n + h n ^ 1 {t n ,y n ; h n ) , with y Q = y , 

where (p 7 denotes some numerical algorithm to compute (p 7 , is called an LLRK scheme. 

When implementing the LLRK discretization (|32j) . that is, when a LLRK scheme is constructed, the 
required evaluations of the expression y„ + cf> (t n ,y n ; .) at t n +\ — t n and c$ (t n +\ — t n ) may be computed by 



different algorithms. In |16j |. 42[ a number of them were reviewed, which yield the following two basic kinds 
of LLRK schemes: 

yn+i =y n + 4> {tn, y n ; h n ) + p(t n , y n ; h n ) , 

and 

y n +\ =z{t n + h n ;t n ,y n ) +/?(t n ,y n ;/i n ) , 
where cj> is a numerical implementation of <f>, z is a numerical solution of the linear ODE 

= B„z(t) + b„ (t) , t e [t n ,t n+ i], (33) 
z(* n ) = yn, (34) 
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and p is the map of the Runge-Kutta scheme applied to the ODE 

dv(t) 



dt 

v (t n ) = 



q(i„,z(t n ) ;t,v(t)), te [t n ,t n+ i], 



(35) 
(36) 



with vector held 



q(in,y«; s,£) = f(s,y„ + <f>(t n ,y n ; s -t n ) + £) - f x (t n ,y n )<f>(t n ,y n ; s - t n ) 
- ft (tn, y«) (s - *n) - f (t n , y n ) > 

for the hrst kind of LLRK scheme, or 

q(in,y n ; = f(s,z(s;i„,y n ) + £) - f x (tn,y n )(z (s;t n ,y n ) - y„) - f t (t„,y„) (s - i„) 

- f (*n,yn) 

for the second one. In the equation (|33p . B„ = f x (t„, y„) is a dxd constant matrix and h n (t) = f t (t n , y n ) (t— 
f n ) + f (t„, y„) — B„y„ is a d-dimensional linear vector function. 

Obviously, a LLRK scheme will preserve the order 7 of the underlaying LLRK discretization only if (f> is 
a suitable approximation to <f). This requirement is considered in the next theorem. 

Theorem 15. Let x be the solution of the ODE with vector field f satisfying the condition H9\) . With 

t n , t n+ i e (t) h , let z n+ i = z„ + h n Ai (t n , z n ; /i„) and v n+ i = v„ + h n A^ (t n , v„; /i„) &e one-step explicit 
integrators of the ODEs t33\) - (34\ ) and i35\) - ![3b)) . respectively. Suppose that these integrators have order of 
convergence r and p, respectively. Further, assume that Ai and A^f fulfill the local Lipschitz condition U6\) . 
Then, for h small enough, the numerical scheme 



y n +i = fn + h n Ai (t n ,y n ; h n ) + h n A% n (t n , 0; h n ) 



satisfies that 



||x(* n+1 )-y n+1 || KChT*™ 
for all tn+i S (t)h> where C is a positive constant. 

Proof. Let X = {x (t) : t £ [to, T]} . Since A" is a compact set contained in the open set DcK'', there exists 
e > such that the compact set 



£ e R d : min |j£ -x(t)|| < e 



is contained in T>. 

First, set y„ = x(t n ) in the equations (f3"3")) - (jM)) and ([33]) - ([56]) . Since x (t n+ i) =yll (t n + h n ; t n ,x(t n )) - 
r (t n + h n ', t n , x(i„)), it is obtained that 



x(i„+i) -x(t„) - h n Ai(t n ,x.(t n );h n ) 
-h n Af tn) (tn,0;h n ) 



< \\<f>(t n ,x(t n ); h n ) - h n Ai(t n ,x(t n ); h n ) 
;t n ,x(t n )) - v(t„+i)[| 



(37) 



where v(i„ + i) is the solution of equation (j35|) - (l36|) at t = 

By definition, r (t n + h n ;t n , x(£„)) is solution of the differential equation 



du (t) 



dt 

u(i„) =0, 



q(t„,x(i„);t,u(i) ), fe [t„,f 
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evaluated at t — t n+ i. Thus, by applying the "fundamental lemma" (see, e.g., Theorem 10.2 in 23[), it is 
obtained that 

||r(i;t n ,x(i n ))-v(i)|| < l(e p C*-*»)_i) (38) 

for t G [t n ,t n+ i], where 

e= sup ||q(t„, x(t„);t, u(t) ) - q(t n ,x(t n );t, u(t) )|| 

tE[t n ,t n+1 ] 

< M \\<p(t n ,x.(t n ); h n ) - h n Ax(t n ,x.(t n ); h n )\\ , 

M = 2 sup ||f x (t, £)||, and P is the Lipschitz constant of the function q(i„,x(i„); •) (which exists by 

Lemma |6]) . 

Furthermore, 



\\(f>(t n ,x.(t n );h n ) - /i„Ai(i n ,x(t n );/i„)|| = ||z(t n+1 ) -z,(t n ) - h n Ax(t n , z (t„) ; h v 



(39) 



since z (i n +i) = x(t„) + (i„, x(i„); /i„) is the solution (|33l) - (|34|) with y„ = x(t„) at t = i n +i- On the other 
hand, 

||z(t n+1 )-z(* n )-/i n Ai(i n ,z(t„);/j n )|| < Cl h r+1 (40) 



and 



v(t n+ i) - /i„A2 ( * n) (tn,0;/i„) 



(41) 



hold, since z n +i — z n + h n K\ (t n , z ra ; /i n ) and v Jt +i — v n + h n b3p (t n , v Tl ; /i n ) are order r and p integrators, 
respectively. Here, c± and c 2 are positive constants independent of h. 
From the inequalities ([37l) - (j4"Tj) . the one-step integrator 



Yn+i = y« + /i n Ai (t n ,y n ; h n ) + h n K\ n (t n , 0; h„) 



has local truncation error 



x(t n +i) -x(i„) - /i n Ai(t n ,x(t n );/i n ) - h n k* {tn \t n) 0]h n ) 



where c = ci + c 2 + ciAf (e p — 1)/P is a positive constant. In addition, since Ai+ A£ with fixed t n , /i n 
is a local Lipschitz function on 2?, Lemma 2 in [|52[ (pp. 92) implies that Ai+ A^*"^ is a Lipschitz function 
on A e C V. Thus, the stated estimate |[x(* n+1 ) -y„+i|| < C/i min{r ' p} for the global error of the LLRK 
scheme y n +i straightforwardly follows from Theorem 3.6 in [23j, where C is a positive contant. Finally, in 
order to guarantee that y n +i £ A e for all n = 0, N — 1, and so that the LLRK scheme is well-defined, it 
is sufficient that < h < S, where S is chosen in such a way that C5 min t r >rf < e. □ 

As an example, consider the computation of the function <j) through a Pade approximation combined 
with the "scaling and squaring" strategy for exponential matrices [2l|. To do so, note that tj> can be written 

as m, 

4> (tn, y«; hn) = Le D "''"r, 

where 

fx(in,y«) ft(t n ,yn) f(tn,y n ) 
1 





D, 



s(d+2)x(d+2) 



L = [ Id 0dx2 ] and r T = [ Qix(d+i) 1 ] i n case °f non-autonomous ODEs; and 



D, 



fx(yn) f(y r , 
o o 



(d+l)x(d+l) 



L = [ Id Odxi 1 a nd r T = [ Oixd 1 ] for autonomous equations. 
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Proposition 16. Set <j> (t n , y n ; h n ) = L (P p . q (2 K ™D„/i„)) " r, where P p , 9 (2 K "D„/i„) is i/ie (p,q)-Pade 



approximation of e 2 



'D„/i„ 



i/ie smallest integer number such that 



< i, and the matrices 



D„,L, r are defined as above. Further, let p be the numerical solution of the ODE h35\) - l3Bjl given by an 
order 7 explicit Runge-Kutta scheme. Then, under the assumptions of Theorem [?| the global error of the 
LLRK scheme 

fn+i = y n + 4>(t n ,y n - h n ) +p(t n ,y n ;h n ) (42) 
for the integration of the ODE fiP-fiP is given by 

||x(*„)-y„|| KMW^i^ 

for all t n £ (t) h , where M is a positive constant. 

Proof. Let JC C T) be a compact set. Since P p . q is an analytical function on the unit circle, it is also a 



Lipschitz function on this region. This and condition 
exists a positive constant L such that 



4 XJjiILji 



< ~ for all t n S (t) h imply that there 



<j>(t n ,&,hn) - <j>(t n ,€i;hn) < ^ 116 

for all £1,62 S /C and i„ £ (t)i. On the other hand, Lemma 4.1 in 43] implies that there exists a positive 
constant M such that 



Z(*n+l) - z {tn) - <t> (tn, 2(t n );h n ) 



< Mh p+q+1 



for all t n € (t) h , where z is the solution of the linear ODE (|3"3")) - (|3~4"1) . 

In addition, since p is an order 7 approximation to the solution of (|3"5|) -([3"o ]) that satisfies the condition 
(ffijl) . the hypotheses of Theorem [T5l hold, which completes the proof. □ 

The next theorem presents a way to define a class of A-stable LLRK schemes on the basis of Pade 
approximations to matrix exponentials. 



Theorem 17. LLRK schemes of the form J^i^ are A-stable if the (p, q)-Pade approximation is taken with 
P < q < P + 2. Moreover, if q = p + 1 or q = p + 2, then such LLRK schemes are also L-stable. 



Proof. Consider the scalar test equation 



dx (t) — Xx (t) dt, 



where A is a complex number with non-positive real part. 

An LLRK scheme of the form (|42l) applied to this autonomous equation results in the recurrence 

+ 1 =y n + 4> (tn, Vn\ h n ) 

= J 7„+L(P p , (? (M)) 2K "r, 



(43) 



where M = 2 Kn T) n h n and 



Here. 



D, 



A \y n 




Pp, 9 (*) = 

denotes the (p, q)— Pade approximation to e z , where 



N P ,,(z) 
D p,-?(^) 



N p>q (z) = 1 + 



V 



-z + 



p(p - 1) 



q + P (q+p)(q+p-l)2\ 
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p(p- l)...l z p 
(q+p)...(q + l) if' 



and D ft? (z) = N g , p (-z). 
Since 



(M) J 



(2^h n xy 





(2- K "/i„A) J y„ 




it can be shown that 



Likewise, 



N P ,,(M) = 



D P ,,(M) = 



N p , 9 (2- K "/i„A) (N M (2-«»/i n A) - 1) y n 



D p , g (2- K "/t„A) (D P) , (2- K -h n X) - 1) y„ 



Hence, 



Therefore, 



Dp~j(M) 



(Dp, g (2- K -/i„A))" 




P P ,,(M) = N P ,,(M)D-;(M) 

N p .,(2-»"h„A) 
D p , g (2-""ft„A) 





(D M (2-^ M )-l) 
(D Pl9 (2-"«h n A)) ^« 
1 



N p , g (2-"^„A) 
D p ,,(2--"/ l „A) 1 I f« 

1 



and so 



(p M (M)r 



N p , g (2~""/t„A) 
D p ,,(2-"»/i„A) 





N p ,,(2~""/i„A) 
D p , g (2-'=«h„A) 



2 K ™ 



- 1 \Vn 



1 



By substituting the above expression in (|43p it is obtained that 

y„+i = R(X)y n , 

where 

i?(A) 



Np,,(2- K -/ l „A) 



D P:? (2-«»/i n A) 

Since 5R(2- K "/i„A) < 0, Theorem 353A, pp. 238 in 6] implies that \R(X)\ < 1 for p < q < p + 2. That is, for 
these values of p and q the LLRK scheme (|42|) is A-stable. The proof concludes by noting that, for q = p+ 1 
or q = p + 2, R{z) = when z — > oo. □ 

From an implementation viewpoint, further simplihcations for LLRK schemes can be achieved in order 
to reduce the computational budget of the algorithms. For instance, if all the Runge Kutta coefficients Cj 
have a minimum common multiple ft, then the LLRK scheme (|42[) can be implemented in terms of a few 
powers of the same matrix exponential e Kh " D " . To illustrate this, let us consider the so called four order 
classical Runge-Kutta scheme (see, e.g., pp. 180 in |6|) with coefficients c = [ , \ 1 ]• This yields 
the following efficient order 4 LLRK scheme 



y«+i = y« + 4>(t ni y n ; h n ) + p(t n ,y n ; h n ) . 



(44) 



where 



P \pn i yVi i hr. 



(2k 2 + 2k 3 + k 4 ), 
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with 



ki = f (t n + Cih n ,y n + cj)(t n ,y n ;Cih n ) + CihrM-x) - f (t„,y„) 

- f x (t n ,y n ) 4> {t n ,y n ; c%h n ) - f* {t ni y n ) Cih n , 

ki = 0, 4>{t n , y„; 4?-) — LAr, 0(t„, y n ; ft„) = LA 2 r, A = (P p , 9 (2- K "D„/i„)) 2K " , 



D, 



fx (^n j Yn ) ft \pn •> Yn ) f(t n ,y n ) 

1 





»(<2+2)x(d+2) 



and K n is the smallest integer number such that 2 Kn D n h n < |. 

Note that the dynamical properties of an order 7 LLRK discretization, as stated in section^ are inherited 
by its numerical implementations if the approximation to the map <fi+ p is o(/i 7-1 ) and smooth enough (i.e., 
of class C 7 ). In particular, these conditions are satisfied by the implementations just introduced, namely, 
those given by (|4"2"|) . This provides theoretical support to the simulation study presented in [3], [l5j], which 
reports satisfactory dynamical behavior of LLRK schemes in the neighborhood of invariant sets of ODEs. 

Finally note that, as an example, this section has focused on a specify kind of LLRK scheme, namely, 
the A-stable scheme (HU) that combines the A-stable Pade algorithm to compute the <p 1 with the 4 order 
classical Runge-Kutta scheme to compute the solution of the auxiliary equation ([55)1 - (1551) . However, because 
of the flexibility in the numerical implementation of the LLRK methods, specific schemes can be designed 
for certain classes of ODEs, i.e., LLRK schemes based on L-stable Pade algorithm and Rosenbrock schemes 
for stiff equations; or LLRK schemes based on Krylov algorithm in case of high dimensional ODEs, etc. For 
all of them the results of this section also apply. 



6. Numerical simulations 

In this section, the performance of the LLRK4 scheme ([4"4"]l is illustrated by means of numerical simula- 
tions. To do so, a variety of ODEs were selected. All simulations were carried out in Matlab2007b, and the 
Matlab function "expm" was used in all computations involving exponential matrices. 

The first example is taken from [2} to illustrate the dynamical behavior of the LLRK4 scheme in the 
neighborhood of hyperbolic stationary points. For comparative purposes, the order 2 Local Linearization 



scheme of 41|, and a straightforward non-adaptive implementation of the order 5 Runge-Kutta formula of 
Dormand & Prince [l8[ (used in Matlab2007b) are considered too. They will be denoted by LL2 and RK45, 
respectively. 
Example 1 

- -2x a + x 2 + 1 - nf (xi , A) , (45) 
^=x 1 -2x 2 + l- M /(x 2 ,A) ! (46) 

where f (u,X) = u (l + u + Xu 2 ) 1 . 

For fi = 15, A = 57, this system has two stable stationary points and one unstable stationary point in 
the region < x±,X2 < 1. There is a nontrivial stable manifold for the unstable point which separates the 
basins of attraction for the two stable points. 

Figure la) presents the phase portrait obtained by the LLRK4 scheme with a very small step-size 
(h = 2~ 13 ), which can be regarded as the exact solution for comparative purposes. The stable manifold M s 
of the unstable point was found by bisection. Figures lb), lc) and Id) show the phase portraits obtained, 
respectively, by the LL2, the RK45 and the LLRK4 schemes with step-size h = 2~ 2 fixed. It can be observed 
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that the RK45 discretization fails to reproduce correctly the phase portrait of the underlying system near 
one of the point attractors. On the contrary, the exact phase portrait is adequately approximated near both 
point attractors by the LL2 and LLRK4 schemes, being the latter much more accurate. Other significant 
difference in the integration of this equation appears near to the stable manifold M s , Changes in the 
intersection point (0,6i) of the approximate stable manifold Mj} with the a^-axis is shown in Table I for 
the considered schemes. The values of 61 were calculated by a bisection method and the estimated order of 
convergence was calculated as 

r — ^ ln( ^ h ~ ^ 2 ) 
hi 2 61/2-61/4 

For h < 2~ 4 , the reported values of r/j for the schemes LL2 and LLRK4 are in concordance with the expected 
asymptotic behavior 61 = ^o + Ch r + 0(h r+1 ) stated by TheoremfTTI and Theorems 3 in [4ljj , respectively, but 
not with the stated by Theorem 3.1 in [2j for the RK45, i.e., w 5. This means that the LL2 and LLRK4 
schemes provide better approximations to the stable and unstable manifolds on bigger neighborhoods of 
the equilibrium points, which is obviously a favorable result for them. These results show out too that the 
LLRK4 scheme preserves much better the basins of attraction of the ODE (|4l)]) - (j4"6]) than the RK45 and 
LL2 schemes. 




Figure 1: Phase portrait of the system H45H - H46I I computed with fixed step-size h by a) LLRK4 scheme with h = 10 -13 (dashed 
line). The unstable point is pointed out with "o"; b) LL2 scheme, with h = 2~ 2 ; c) RK45 scheme, with h = 2 -2 ; d) LLRK4 
scheme, with h = 2 -2 . In all cases the solid lines represent the solution computed with h = 2~ 2 . 



In what follows, we compare the accuracy of the LLRK4 scheme with those of the LL2 scheme, and 
the Matlab2007b codes ode45 and odel5s in the integration of a variety of ODEs. We recall that the code 
ode45 is a variable step-size implementation of the explicit Runge-Kutta (4, 5) pair of Dormand & Prince 
[l8j ]. which is considered for many authors the most recommendable scheme to apply as a first try for most 
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problems. On the other hand, the code odel5s is a quasi-constant step-size implementation in terms of 
backward differences of the Klopfenstein-Shampine family of numerical differentiation formulas of orders 
1 — 5, which is designed for stiff problems when the ode45 fails to provide desired result [551 ] . 



Step-size\Scheme 


LL2 


RK45 


LLRK4 


h 


a 


rh 


& 


rh 




rh 




0.71911 




0.70377 




0.56615 




2- 2 


0.69688 


1.931 


0.53673 


3.00 


0.58441 


5.142 


2- 3 


0.61727 


2.190 


0.59859 


4.18 


0.59032 


2.384 


2- 4 


0.59639 


2.145 


0.59088 


7.23 


0.59049 


3.354 


2- 5 


0.59182 


2.056 


0.590458 


6.93 


0.590459 


3.901 


2- 6 


0.59079 


2.027 


0.59045593 


6.39 


0.5904561 


3.973 


2- 7 


0.59054 


2.014 


0.5904559168 


5.98 


0.590455917 


3.989 


2 -8 


0.59048 




0.5904559165 




0.590455916 





Table I. Values of and rh computed by the LL2, RK45 and LLRK4 schemes in the integration of the system 
[|3g ]> -([36 |> . for different values of h. 



In order to compare the (non-adaptive) LL schemes with the adaptive Matlab codes, the following 
procedure was carried out. First, one of the Matlab codes is used to compute the solution with fixed values 
of relative (RT) and absolute (AT) tolerance. Then, the resulting integration steps (t)h are set as input in 
the other schemes for obtaining solutions at the same integration steps. Second, the Matlab code odel5s 
is used to compute on (t)h a very accurate solution z with RT — RA — 10~ 13 . Third, the approximate 
solution y of the ODE is computed for each scheme on (t)h, and the relative error 



RE 



max 

--l,...,d; tjG(t) h 



i(tj)-yi(tj) 



is evaluated. 

The following four examples are of the form 



dt 



= Ax + f(x), 



(47) 



where A is a square constant matrix, and f is a nonlinear function of x. The vector field of the first two ones 
has Jacobians with eigenvalues on or near to the ima gina ry axis, which make these oscillators difficulty to 
be integrated by a number of conventional integrators [20j, 1 591 . The other two are also hard for conventional 
explicit schemes since they are examples of stiff equations }55| . Example 5 has an additional complexity for 
a number of integrators that do not update the Jacobians of the vector field at each integration step [H^, H^] : 
the Jacobian of the linear term has positive eigenvalues, which results a problem for the integration in a 
neighborhood of the stable equilibrium point x = 1. 
Example 2. Periodic linear: 



rfx 



- = A(x + 2) 



with 



i 
-i 



xi(t ) - -2.5, x 2 (i ) = -1.5, and [t ,T] = [0,4tt]. 
Example 3. Periodic linear plus nonlinear part: 



rfx 
~dt 



A(x + 2) 



O.lx 2 , 
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where the matrix A is defined as in the previous example, x(i ) = 1, and [to,T] = [0,47r]. 
Example 4. Stiff equation: 

| = -100H(x + l), 

where H is the 12-dimcnsional Hilbcrt matrix (with conditioned number 1.69 x 10 16 ), Xj(to) = 1, and 
[t ,T] = [0,1]. 

Example 5. Stiff linear plus nonlinear part: 

d~x 

^ = 100H(x - 1) + 100(x - l) 2 - 60(x 3 - 1), 
where H is the 12-dimensional Hilbcrt matrix, x$(io) = —0.5, and [to,T] = [0, 1]. 



Example 


Scheme 


Relative 
Tolerance 


Absolute 
Tolerance 


TVTC! 

JNo 


Relative 
Error 




odel5s* 


io- 3 


io- b 


334 


0.19 


2 : Periodic linear 


odeAb 
LL2 


5 x 10~ 6 


5 x 10~ 9 


340 
334 


8.2 x 10" 5 
1.6 x 10~ 12 




LLRKA 






334 


1.6 x 10~ 12 




odel5s* 


nr 3 


10" e 


287 


0.30 


3 : Periodic linear 


odeAb 


6 x 10~ 6 


6 x 10~ 9 


289 


1.3 x 10~ 4 


plus nonlinear part 


LL2 






287 


3.1 x 10~ 2 




LLRKA 






287 


1.1 x 10" 5 




odel5s* 


HT 3 


IO"* 3 


66 


6.7 x 10" 2 


4 : Stiff linear 


odeAb 
LL2 


5 x 10~ 4 


5 x 10~ 7 


66 
66 


5.3 x 10" 3 
1.8 x 10~ 10 




LLRKA 






66 


1.8 x 10" 10 




odelbs* 


10~ 2 


10" 4 


49 


0.31 


5 : Stiff linear 


odeA5 


nr 1 


10~ 3 


104 


0.37 


plus nonlinear part 


LL2 






49 


0.43 




LLRKA 






49 


4.3 x 10~ 5 




odel5s 


10~ 2 


10" b 


103 


0.35 


6 : Nonlinear 


odeA5* 


10~ 3 


10" 6 


47 


0.08 


(no stiff) 


LL2 






47 


4.19 




LLRKA 






47 


0.25 




odel5s 


1.5 x 10 _a 


1.5 x wr ri 


2281 


1.2 x 10~ 3 


7 : Nonlinear 


odeAh* 


io- 7 


io- 10 


2285 


1.6 x 10~ 3 


(moderate stiff) 


LL2 






2285 


400 




LLRKA 






2285 


6.9 x 10~ 3 



Table II. Accuracy of the LL2, LLRK4, ode45 and odel5s schemes in the integration of examples (2)-(7). With 
the symbol * is denoted the Matlab code used to set the time partition (t) h in each example. NS denotes the number 
of steps required for each scheme to compute the solution on (t) h . 

The results of the integration of these equations for each scheme are shown in Table II. For illustration, 
Figure 2 shows the path of the variable Xi and its approximation yi obtained by the LLRK4 scheme in the 
integration of these equations. Remarkable, in all the examples, the relative error of the solution obtained by 
the LLRK4 scheme is much lower that those of the LL2, ode45 and odel5s with the same or lower number of 
steps. These results are easily comprehensible for five reasons: 1) the dynamics of these equations strongly 
depend on the linear part of their vector fields; 2) the LL2 and LLRK4 schemes preserve the stability of 
the linear systems for all step-sizes, which is not so for conventional explicit integrators; 3) the LL2 and 
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LLRK4 schemes are able to " exactly" (up to the precision of the floating-point arithmetic) integrate linear 
ODEs, which is a property not satisfied by for conventional explicit and implicit schemes; 4) the LL2 and 
LLRK4 schemes update the exact Jacobian of the vector field at each integration step, which is not done 
by most of conventional schemes; and 5) the LLRK4 has higher order of convergence than the LL2 scheme. 
Further, note that although the LLRK4 scheme is not designed for the integration of stiff ODEs in general 
(because the auxiliary equation (l3"5jl - (131)1) might "inherit" the stiffness of the original one) it is clear that, 
by construction, it is suitable for equations with stiffness confined to the linear part. Example are the 
classes of stiff linear and semilinear equations represented in the Examples 2 and 3. This is so, because at 
each integration step the stiff linear term is locally removed from the vector field of the auxiliary equation 
(|35|) and, in this way, the stiff linear part is well integrated by the (A-stable) LL scheme and the resulting 
non-stiff equation (l3"5j) can be well integrated by the explicit RK scheme. 

The following two examples are well known nonlinear oscillators. 

Example 6. Non-stiff nonlinear: 



rfxi 
~dT 



= 1 + x?x 2 - 4xi, 



d*2 „ 2 

— = 3x x - x lX2 

where xi(to) = 1.5, x 2 (io) = 3, and [to,T] = [0,20]. This equation, known as Brusselator equation, 
typical test equation of non-stiff nonlinear problems (see, e.g., 23] ) . 
Example 7. Mild-stiff nonlinear: 



dxi 

= e((l - x2)xi + x 2 ) 



dt X2 ' 
dt 



where e = 10 , xi(£o) = 2, X2(io) = 0, and [to,? 1 ] = [0,2]. This equation, known as Van der Pol equation, 



is a typical test equation of stiff nonlinear problems (see, e.g., [24| ). 

The results of the integration of last two equations for each scheme are also shown in Table II and Figure 
2. For these equations, the relative error of the solutions obtained by the LLRK4 scheme is much lower 
that those of the LL2, but quite similar to those of the codes ode45 and odel5s (which have higher order 
of convergence). This indicates that the LLRK4 scheme is also appropriate for integrating non-stiff and 
mild-stiff nonlinear problems as well. 

In summary, results of Table II clearly indicate that the non adaptive implementation of the LLRK4 
scheme provides similar or much better accuracy than the Matlab codes with equal or lower number of 
steps in the integration of variety of equations. This suggests that adaptive implementations the LLRK 
discretizations might archive similar accuracy than the Matlab codes with lower or much lower number of 
steps, a subject that has been already studied in 5|| 44 1. 



Finally, we want to point out that equations of type (|47p frequently arises from the discretization of 
nonlinear partial differential equations. In such a case, mild or high dimensional ODEs of that form are 
obtained and, as it is obvious, LLRK schemes like ([4"4"]l based on Pade approximations are not appropriate. 
Nevertheless, because the flexibility of the high order Local Linearization approach described in Section [2j 
feasible high order LL schemes can be designed for this purpose too. For instance, by taking into account 
that ^ 

(t>(t n ,y n ; y) = <^(yfx(y«))f(yn), 

where (p(z) — (e z — l)/z, the LLRK4 scheme (|44| can easily modified to defined an order 4 LLRK scheme 
for high dimensional ODEs. Indeed, such scheme can be defined by the same expression ((44]), but replacing 
the formulas of <p(t n ,y n ; 4f ) and 4>{t n ,y n ; h n ) by 

<j>(tn,y n ; y) = ^(yfx(y n ))f(y«) 
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example 2 



example 3 






10 12 14 




0.5 



example 5 



0.2 0.4 0.6 0.8 1 

t 

example 7 



1.5 



Figure 2: Path of the variables xi (solid line) and its approximation yi (dots) obtained by the LLRK4 scheme in the integration 
of the ODEs of examples 2-7. The time partition (i)^ used in each case for yi is specified in Table II. The "exact" path of xi 
is computed with the Matlab code ode!5s with RT = RA = 10~ 13 on a very thin partition. 
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and 

<j>(t n ,y n ;h n ) = ^^-fx(yn)^(Y fx (y»)) + I ) #™>y«; 

resp ectively, where tp denotes the approximation to tp provided by the Krylov subspace method (see, i.e., 
[271]). Then, a comparison with exponential- type integrators designed for high dimensional equations of the 
form (|47p can be carried out, but this subject is out of the scope of this paper. 



7. Conclusions 

In summary, this paper has shown the following: 1) the LLRK approach defines a general class of high 
order A-stable explicit integrators; 2) in contrast with others A-stable explicit methods (such as Rosenbrock 
or the Exponential integrators) , the RK coefficients involved in the LLRK integrators are not constrained by 
any stability condition and they just need to satisfy the usual, well-known order conditions of RK schemes, 
which makes the LLRK approach more flexible and simple; 3) LLRK integrators have a number of convenient 
dynamical properties as the linearization preserving and the conservation of the exact solution dynamics 
around hyperbolic equilibrium points and periodic orbits; 4) unlike the majority of the previous published 
works on exponential integrators, the above mentioned convergence, stability and dynamical properties are 
studied not only for the discretizations but also for the numerical schemes that implement them in practice; 
5) because of the flexibility in the numerical implementation of the LLRK methods, specific-purpose schemes 
can be designed for certain classes of ODEs, e.g., for stiff equations, high dimensional systems of equations, 
etc.; 6) order 4 LLRK formula considered in this paper provides similar or much better accuracy than the 
order 5 Matlab codes with equal or lower number of steps in the integration of variety of equations, as well 
as, much better reproduction of the dynamics of the underlying equation near stationary hyperbolic points. 

Finally, it is worth to point out that theoretical properties of the LLRK methods studied here stro ngly 
support the results of the numerical experiments carried out by the authors in previous works 14 1, [15| . 



in which the performance of other LLRK schemes is compared with that of existing explicit and implicit 
schemes. 
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